function [H,P,Hmax,backP] = sim_supply(Rscale,q0,Theta,cost)
run set_parameters.m

% for age-dependent polynomial h(a), degree K = 4
X_age = [ones(T,1) age' (age.^2)' (age.^3)' (age.^4)'];
h_age = X_age*Theta;
h_age(11:end) = h_age(10);

sum_ha = zeros(T,1); % sum of h(a) from age 1:a
for a = 1:T
    sum_ha(a) = (age <= a)*h_age;
end


% Parameters in Recruitment function
Rmax = 5277994; % Maximum BFT Recruitment in sample
k = Rscale*Rmax;
rho = 8.1;

Emax = 10000; % set maximum effort
Elen = 10000; % length of the effort vector
E = linspace(0,Emax,Elen); % effort (range from 0 to Emax)

% record simulated values
N = zeros(Elen,T); % population in each age class
N0 = zeros(Elen,1); % recruitment
H = zeros(Elen,1); % harvest
P = zeros(Elen,1); % price

for i = 1:Elen % for each effort level
    e = E(i);
    % recruitment N0: function of e
    Nb = (mature'.*w_a)*exp(-q0*e*sum_ha - sum_m);
    N0(i) = k - (k/(rho*Nb));
    % population N(a): function of N0 & e
    N(i,:) = N0(i)*exp(-q0*e*sum_ha' - sum_m');
    % harvest H: function of e -------------------------
    H(i) = N(i,:)*(q0*e*h_age./(q0*e*h_age + m_a).*(1-exp(-q0*e*h_age-m_a)).*w_a');
    % price P: function of e ---------------------------
    P(i) = cost*e./H(i);
end % for i

% get the backward bending Price & Harvest
[Hmax, maxIndx] = max(H);
backP = P(maxIndx);
end